Hydrogen Electron Capture in Accreting Neutron Stars and the Resulting 

g-Mode Oscillation Spectrum 



Lars Bildsten and Andrew Gumming 
Department of Physics and Department of Astronomy 

University of California, Berkeley, CA 94720 
bildsten@fire.berkeley.edu, cumming@fire.berkeley.edu 



ABSTRACT 

We investigate hydrogen electron capture in the ocean of neutron stars accreting at rates 
10^^" <^ M <^ 10~^ Mq yr~^. These stars burn the accreted hydrogen and helium unstably 
in the upper atmosphere {p < 10^ g cm~^) and accumulate material which usually contains 
some small amounts of hydrogen (mass fractions are typically Xr ~ 0.1 — 0.2) mixed in with the 
heavier iron group ashes. The subsequent evolution of this matter is determined by compression 
towards higher densities until electron capture on the hydrogen occurs. We construct steady-state 
models of the electron captures and the subsequent neutron recombinations onto the heavy nuclei. 
The density discontinuity from these captures gives rise to a new g-mode (much like a surface 
wave), which has a lowest order {I = 1) frequency of ~ 35 Hz(Xr/0.1)^/^ on a slowly rotating 
{fs ^ 30 Hz) star. We also discuss, for the first time, a new set of non-radial g-modes unique to 
these high accretion rate neutron stars. These modes are "trapped" in the finite thickness layer 
where the electron captures are occurring. The lowest order mode frequencies are in the 1 — 10 Hz 
range for a few radial nodes on a slowly rotating star. Though the majority of the mode energy 
resides in the electron capture transition layer, the eigenmode propagates to higher altitudes 
above the layer and can thus be potentially observable and excited by the nuclear burning or 
other mechanisms. We also show that the density jump splits the ocean's thermal modes into 
two distinct sets, which have most of their energy either above or below the discontinuity. We 
conclude by discussing how the dispersion relations for these modes are modified for a rapidly 
rotating {fs ^ 30 Hz) neutron star. Whether any of these modes are observable depends on the 
damping mechanism and the ability to excite them, issues we will address in a future paper. 
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1. Introduction 

The liquid core, solid crust and overlying ocean and atmosphere give a neutron star (NS) a very rich 
spectrum of non-radial oscillations, with acoustic p- modes in the 10 kHz range and a variety of lower frequency 
g-modes from the toroidal and spheroidal displacements of the crust and surface layers (sec McDermott, Van 
Horn & Hansen 1988 for an overview). Non-radial g-mode oscillations have been studied extensively for 
isolated and slowly rotating radio pulsars, where the atmosphere and ocean are quiescent. These studies can 
be separated by the sources of buoyancy that were considered: (1) entropy gradients (McDermott, Van Horn 
& Scholl 1983; McDermott et al. 1988), (2) density discontinuities (Finn 1987; McDermott 1990; Strohmayer 
1993) and (3) mean molecular weight gradients due to /3-equilibrium in the core (Reisenegger & Goldreich 
1992). Unfortunately, none of these modes have been securely identified with any particular radio pulsar 
phenomenon. 

There arc also many accreting neutron stars in binary systems that emit much of their accretion lumi- 
nosity in X-rays (Lewin, van Paradijs & van den Heuvel 1995). The thermal and compositional structure 
of these accreting NSs is quite different from that of the isolated radio pulsars, and is far from equilibrium. 
There have been previous studies of non-radial oscillations in these objects (McDermott & Taam 1987; Bild- 
sten & Cutler 1995 (hereafter BC95); Bildstcn, Ushomirsky & Cutler 1996 (hereafter BUC96); Strohmayer 
& Lee 1996 (hereafter (SL96)). In this paper, we calculate the g-mode spectrum of the ocean of an accreting 
NS, focusing on the critical effects of hydrogen electron capture. 

1.1. Observational Motivation 

Much of this work is motivated by the rich spectrum of quasi-periodic oscillations (QPOs) seen in the 
luminosity of the brightest accreting neutron stars. Utilizing EXOSAT data., Hasinger & van der Klis (1989) 
found that the highest accretion rate objects {M > 1O~^°M0 yr~^) fall into two separate classes. Six 
objects trace out all or part of a "Z" in an X-ray color-color diagram and exhibit time-dependent behavior 
that correlates with the position along the Z. These "Z sources" have QPO's in the 15-50 Hz and 5-7 Hz 
range with up to 10% modulation of the X-ray flux. The other objects fall into separated regions of the 
color-color diagram and do not show this distinctive QPO phenomenology. These "Atoll sources" accrete at 
lower rates than the Z sources and exhibit Type I X-ray bursts resulting from the unstable ignition of the 
accumulated hydrogen and helium (see Bildsten 1998a for a recent review). 

The Rossi X-Ray Timing Explorer (RXTE) has recently added to this list of QPO phenomenology. 
Two drifting kHz QPOs are seen in the persistent emission from both Atoll and Z sources, and coherent 
oscillations have been observed during Type I X-ray bursts from several Atoll sources (see Strohmayer et 
al. 1996 for an example). In these objects, the kHz QPOs are separated by a frequency identical to that 
seen in the burst, leading naturally to a beat frequency model in which the difference frequency is the NS 
spin frequency (Strohmayer et al. 1996). The physical origin of the upper frequency differs in the models 
(see van der Klis 1998 for a summary), and usually involves the Kepler frequency at some special place in 
the accretion disk (Kaaret, Ford & Chen 1997; Zhang, Strohmayer & Swank 1997; Miller, Lamb & Psaltis 
1998). The inferred spin frequencies are in a remarkably narrow range, 250-500 Hz (White & Zhang 1997; 
van der Klis 1998; Bildsten 1998b). There are also many QPOs with frequencies less than 100 Hz seen in 
the Atoll sources. Stella & Vietri (1998) proposed that these might be due to Lense-Thirring precession of 
the inner accretion disk. 

It may be that all of the < 100 Hz and kHz QPOs can be accommodated within the existing models 
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(see van dcr Klis 1995, 1998) that use the accretion disk and/or spherical flow to generate the periodic 
phenomena. However, it is important to pursue non-radial oscillations of the surface layers of the NS as a 
possible source of some of these periodicities. The predicted frequencies are well-constrained, depend on the 
underlying NS structure and have well-understood dispersion relations, even when the NS is rapidly rotating 
(BUC96). This means that the successful identification of a non-radial pulsation with an observed frequency 
would tell us much about the thermal and compositional makeup of the NS and measure or constrain the 
rotation rate and magnetic field. 

1.2. Accreting Neutron Star Structure and Non-Radial Oscillations 

In most cases where stellar pulsations are studied, the underlying stellar model is reasonably well 
understood and constrained via other observations. Initial theoretical work then consists of finding the 
adiabatic mode structure and exploring possible excitation mechanisms, usually focusing on the K-mechanism 
which relies on opacity changes during the pulsations (see Cox 1980 for an excellent overview). The situation 
is quite different in an accreting NS. The conditions in the outer layers can change on short timescales, 
the deep internal structure depends on the accretion history and there is a broader range of excitation 
possibilities. We now briefly outline the underlying NS structure and the resulting g-mode spectrum. 

Figure |] is a schematic of the NS atmosphere and ocean. For accretion rates M < IQ^^Mq yr^^, the 
accreted hydrogen (H) and helium (He) rich material burns unstably at a density of p ^ 10^ g cm^^ and 
within a few hours to days upon arrival on the star. The very high temperatures reached (T > 10^ K) 
during the thermal instability produce elements at and beyond the iron group. The isotopic mixture from 
this burning is not well known, though it seems to always be the case that a substantial amount of H (the 
residual mass fraction is typically Xr ~ 0.1) remains unburned. The ashes from the burning accumulate on 
the NS, undergo further compression and form a relativistically degenerate ocean. The hydrogen is eventually 
depleted by electron captures at p ~ 2 x 10 g cm~^, leading to an abrupt rise in density at that location. 
At still higher densities, p ~ 10^-10^° g cm~'^, the material crystallizes and forms the crust. The thickness 
of the ocean is « 10*^ cm ^ R. There is no overt evidence for magnetic fields on these NSs (i.e. none of them 
are persistent X-ray pulsars) and most arguments about the nature of Type I X-ray bursts limit B < 10^ G 
(Bildsten 1998a), weak enough not to affect the ocean g-modes (BC95). 

There has been prior work on g-modes in the upper atmosphere of the NS, above and around the H/He 
burning region. McDermott & Taam (1987) calculated g-modes of a bursting atmosphere. SL96 calculated 
the non-adiabatic mode structure for atmospheres accreting and burning in steady-state and found that g- 
modes may be excited by the e-mechanism (i.e. by the nuclear burning) when M < 10~^'^Mq yr~^. However, 
steady-state burning does not occur at these low M's, and realistic calculations for the time dependent NS 
atmosphere have yet to be carried out. 

Our focus in this paper is on the adiabatic mode structure in the deep ocean underneath the H/He 
burning where the thermal time (hours to days) is much longer than the mode period. The g-modes reside in 
the ocean as they are effectively excluded from the crust by the restoring force from the finite shear modulus 
there (McDermott et al. 1988; BC95). If the NS is slowly rotating, the dispersion relation for these modes 
is that for shallow water waves, oc k^, where k = + l))^/^/i? is the transverse wavenumber and / is 
the angular eigenvalue. The transverse wavenumber changes if the star is rapidly rotating, but this is easily 
calculable (BUC96) and it is straightforward to find the g-mode spectrum from the non-rotating calculations. 

The different sources of buoyancy yield a rich spectrum of g-modes. The abrupt rise in density associated 
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with the hydrogen electron capture boundary layer supports a density discontinuity mode of frequency 
fd ~ 35 Hz(Xr/0.1)^/^ for Z = 1 and R = 10 km (equation The internal buoyancy due to the 

composition gradient within the electron capture boundary layer creates a new spectrum of modes which 
are "trapped" . Much of the mode energy is confined to the boundary layer and our WKB estimate of the 
mode frequency gives ftr ~ 8.5 Hz/ntr(^r/0.1)^^^ for / = 1 and R — 10 km (equation where ritr is 

the number of nodes in the boundary layer. There is also a set of thermal g-modes (BC95) in the same 
frequency range as the trapped modes, which are separated by the density discontinuity into two distinct 
sets (see §4.3). 



1.3. Outline of This Paper 

We begin in §2 by explaining the physics of the deep ocean, the electron captures that occur there 
and the subsequent neutron captures by the heavy nuclei. We introduce a convenient fitting formula for the 
electron capture rate and use it to generate a series of stellar models for the subsequent seismic investigation. 
Section 3 summarizes the equations for non-radial oscillations, discusses the various sources of buoyancy in 
the deep ocean and our method for solving the eigenvalue problem. We present the full spectrum of g-modes 
for isothermal oceans in §4. There, we show both the modes due to the electron capture discontinuity and the 
"trapped" modes, and find that the thermal modes are split into two distinct sets by the density discontinuity. 
We use the WKB approximation to obtain reliable frequency estimates for these modes. We construct a 
few realistic non-isothermal ocean models in §5 and show their g-mode spectra. We summarize our results 
in §6 and discuss the effects of rapid rotation and the important issues of mode excitation, damping and 
transmittal of the signal to the observer. 



2. Electron Capture on Hydrogen in Accreting Neutron Stars 

Rosenbluth et al. (1973) first noted that accreted matter would undergo electron capture on hydrogen 
(e~ + p —> n + Ve) once high enough densities were reached (p > lO'^ gcm~^) for the electron Fermi energy 
to exceed the 1293 keV threshold. Hydrogen electron captures have since been studied in two very different 
contexts in accreting neutron stars. The first is at very low accretion rates, where the thermonuclear 
consumption of hydrogen (say via the CNO or pp cycle) is so slow that the matter reaches the capture 
density prior to having burned any hydrogen. In this hydrogen-rich environment, the neutrons created by 
this process are immediately captured by protons (p + n ^ D -|- 7) . Further reactions burn this matter to 
at least helium, releasing 7 MeV per accreted proton. Van Horn & Hansen (1974) modeled transient X-ray 
sources as the result of electron capture triggered thermonuclear fiashes on low mass (M < 0.15Mq) neutron 
stars. A variant of this model was later invoked as a possible source of local gamma-ray bursts by Hameury et 
al. (1982), who found that the accretion rates where the electron capture ignites a thermonuclear instability 
is comparable to that expected from a lone neutron star steadily accreting from the interstellar medium. 

The second context is at the high accretion rates typical for the brightest X-ray sources. In these stars, 
the matter undergoes unstable H/He burning within a few hours to days of arrival on the NS, long before 
reaching densities high enough for electron capture. Provided M > Mc2 ~ lO^^^M© yr"-'^(2'cNo/0.01)"'^/^ 
(Bildsten 1998a), the H burning (via the /^-limited hot CNO cycle) is thermally stable while the H accumu- 
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lates. For accretion rates in excess of 

Afa » 7 X 10->M,y.- (h^y f-^Y" (^)""\ (1) 

"■^ V ^ / V 10 km J V 0-01 / 

(Bildsten 1998a) the CNO cycle cannot consume all of the H before the He unstably ignites (Lamb & 
Lamb 1978; Taam & Picklum 1979; Fujimoto, Hanawa & Miyaji 1981) in a local thermonuclear flash. For 
Mc2 < M < Mci, the H is completely burned in a stable manner and the accumulated pile of pure He 
unstably ignites. There are many neutron stars accreting at rates comparable to Mci- However, the strong 
dependence of Aid on Zcno smd R makes it difficult to know in which regime a particular star is burning, 
especially if a large amount of spallation occurs when the infalling matter decelerates at the surface (Bildsten, 
Salpeter & Wasserman 1992). [] The ashes from these two types of unstable burning are quite different. In the 
pure He case, the ashes are mostly ^^Ni (Joss 1978; Ayasli & Joss 1982) and there is little to no hydrogen left 
at densities in excess of 10^ g cm^'^. Things are different when the He ignites and burns in a hydrogen- rich 
environment (M > Md), as we now discuss. 

During the thermonuclear flash in the mixed H/He regime, hydrogen burning is accelerated by freshly 
minted seed nuclei produced by He burning, enhancing the local energy release and substantially complicating 
the nuclear reaction chains. The very high temperatures reached (T ;> 10^ K) during the thermal instability 
produce elements beyond the iron group (Hanawa, Sugimoto & Hashimoto 1983; Wallace & Woosley 1984; 
Schatz et al. 1997). The isotopic mixture from this burning is still not well known, though everyone 
agrees that a substantial amount of hydrogen (the residual mass fraction is typically Xr ^ 0.1-0.2) remains 
unburned (Ayash & Joss 1982; Woosley & Weaver 1984; Taam et al. 1993; Schatz et al. 1997). As the 
mix of H and extremely heavy elements accumulates beneath the H/He burning layers, at temperatures 
in the range T « (2-5) x 10* K, some proton captures by the heavy nuclei occur (Schatz 1997, private 
communication). However, the main proton consumption occurs when the matter has been compressed by 
accretion to densities high enough for hydrogen electron capture (Woosley & Weaver 1984; Fushiki et al. 
1992; Taam et al. 1996). Electron capture onto the nuclei themselves occurs at much greater densities and 
has been discussed previously by Sato (1979), Haensel & Zdunik (199Ga, 1990b) and Blaes et al. (1990). 

In this section, we discuss the nuclear physics of the electron capture transition layer in the ocean 
underneath the H/He burning layer for those objects with M > Mci- We find that the neutrons produced 
by the electron captures recombine with the heavy nuclei, substantially changing the mix of heavy elements 
in the layer and below. Our present motivation is to calculate the internal buoyancy associated with the 
composition gradient in the layer where electron captures are occurring, and the resulting spectrum of g- 
modes. We find trapped g-modes and a density discontinuity g-mode that directly depend on the internal 
buoyancy of the electron capture transition layer, and thus only arise when hydrogen is present at these large 
depths. The presence of such g-modes would identify the star as being in the M > Md regime. 



'^Many of the recent interpretations of the kHz QPO's from these same objects imply that the NS radius is inside the last 
stable orbit (Kaaret et al. 1997; Zhang et al. 1997), resurrecting some of the original suggestions of Kluzniak & Wagoner (1985) 
and Kluzniak, Michelson & Wagoner (1990). If true, then spallation might reach the levels found by Bildsten et al. (1992), as 
the flow can have a large radial component when it finally reaches the surface. 
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2.1. The Electron Capture Rate for Hydrogen 

The electron capture rate for hydrogen (i.e. + p ^ n + v^,) residing in an electron gas at temperature 
T with Fermi energy Ep (including electron rest mass) is 

Rec{EF,T)^{^^^I{EF,T), (2) 

(Lang 1980) where ft = 1065 s as inferred from the measured neutron half-life (Barnett et al. 1996) and the 
known / value. The function / is a dimensionless phase space integral over the electron energy i?, 

T(F T)- ^ r E{E^-mlc^Y/\E-QfdE 
^ ^' ' {m.c^f Jq l + eM{E-Ep)/kBT] ' 

where Q = {nin — mp)(? = 1293.318 keV is the energy threshold for the reaction. We have evaluated the 
integral in equation (|^) and agree with the tabulated rates of Oda et al. (1994) to within 2%. For the 
purposes of both analytical insight and numerical speed, we now discuss a few limiting forms for the capture 
rate and present a convenient fitting formula. 

Much insight is gained if we assume relativistic electrons and neglect the m-gC^ term in the integral J. 
There are then two limits to discuss. For T = 0, captures only occur once the threshold is exceeded. For 
values of the Fermi energy slightly above threshold {Ep — Q ^ Q), the integral becomes 

I{Ef.T) - (-^\ / -Ef - Q \ ^ for T = near threshold. (4) 



3 \m!,c^ J \ Q 

When Ep < Q, captures only occur if the electrons have a finite temperature. The integral can be approxi- 
mately found when Q — Ep 3> ksT, giving 

I{Ep, T) — --^ exp , for Q - > ksT, (5) 



(meC^)5 \ ksT 

so that most captured electrons are far out on the thermal tail of the distribution. 

These limiting forms motivate a fitting formulae for the rates that we now introduce. Since this might 
prove a convenient form for others to use, we motivate our choice and discuss how to use it. It represents 
a modified (and smooth) combination of equations (^) and (^. The electron capture rate is temperature 
sensitive even above threshold when Ep — Q ~ ksT. We thus add a term oc fcsT to Ep in equation (^) which 
provides some thermal "fuzziness" to the Fermi energy. The proportionality constant and Fermi energy to 
connect the modified equation (^ with the unmodified equation (|^) are found by demanding that the rate 
and its first derivative with respect to Ep (at constant T) are continuous there. This results in the form 

J ^, [ 2{kBTf exp[{Ep - Q)/kBT] Ep < Q + kBT\n{%/2) 

^ (mec2)5 ^ \ {l/?,){Ep + (3 - In9/2)A:BT ~ Qf Ep > Q + fcBrin(9/2), ^ ' 

which reproduces the electron capture rate within a factor of two over a large range of temperatures and 
densities. Figure ^ shows the ratio of the true integral (eqn. |^ ) to that given by equation (^ for a range 
of Ep and temperatures of Tg = 1, 2, 5, 8 and 10. Clearly the fit does rather well at representing a function 
that is changing by orders of magnitude in this range. The growing deviation at Ep ^ Q arises from our 
neglect of the next order term in equation (||), which is oc {Ep — Q)'*. 

We implement this fit by first tabulating (as a function of Ep and T) the ratio of the real answer to the 
fit. Then, for any Ep and T, we use equation (g) to evaluate / and then find the correction via interpolation 
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within this prc-made table. This proves to be quite accurate with a relatively small table, as the fitting 
formula takes care of most of the strong variations in /. 



We describe the structure of the hydrogen electron capture transition layer in j |2.3| , but we can obtain 
a rough picture of where the electron captures occur by noting that in steady state, the hydrogen depletes 
when the lifetime to electron capture, tec = ^/Rec, becomes comparable to the time to reach that depth, 
taccr = 11/1^, where m is the local accretion rate (mass accreted per unit time per unit area) and y ~ p/g is 
the column depth. The local Eddington accretion rate is 

mEdd = = 7.5 X 10* Me (t7^) g cm-^ s-\ (7) 



arR V 10 ^™ 

where ax is the Thomson cross-section and /ig is the mean molecular weight per electron. For simplicity, 
we take niEdd to be the value for /ig = 1 and R = 10 km. Figure ^ shows tec as a function of Ep for many 
different temperatures. From top to bottom, the solid lines are for Tg = 2, 3, 4, 5, 6, 7, 8, and 9. The vertical 
dashed line displays Ep = Q- The lower (upper) dotted line shows iaccr for rh — m^dd ("^ — O.lmEdd), 
where, for simplicity, we have assumed the degenerate electrons provide all the pressure. It is clear that 
for Tg ^ 6 at m = m^dd and Tg ^ 3 at m = O.lwEdd, substantial pre-threshold capture occurs, i.e. many 



protons are consumed while Ep < Q. We discuss this further in §2.3 



2.2. The Fate of the Free Neutrons 

In the hydrogen- rich environment, the neutrons produced by the electron captures rapidly thcrmalize 
via elastic scattering with protons and undergo a radiative capture (free decay is irrelevant for these high 
densities, as the radiative capture times are <C 10^^° seconds). Woosley & Weaver (1984) and Taam et al. 
(1996) assumed that the neutrons would capture onto protons, starting a chain of rapid nuclear reactions up to 
the Fe group. We find that this is not the case. The reaction rate for capture onto protons (n -I- p — > D 7) 
is {(jv) — 7.3 X 10~^° cm^ s~^ (Caughlan & Fowler 1988), whereas typical rates for neutron captures onto 
heavy nuclei are 10~^^-10~^^ cm"^ s~^, nearly two orders of magnitude larger. We thus find that, for the 
hydrogen mass fractions we consider, most neutrons capture onto the heavy nuclei, releasing the typical 7-8 
MeV binding energy of a nucleon in a heavy nucleus. For example, from Cowan, Thielemann & Truran 
(1991), the neutron capture rate for ^^Fe is 5.8 x lO^^^cm'^ (at 30 keV), so that neutrons will capture 
preferentially onto iron rather than protons when their ratio by number is < 80, requiring X < 0.58. Similar 
conclusions apply to most nuclei we could consider. 

The neutron captures substantially modify the isotopic mix of heavy elements. At a particular depth in 
the transition layer, those elements with larger neutron capture cross sections grow more and more neutron 
rich mitil they become /3-unstable on a timescale comparable to the fiow time across the region (i.e. days, 
see Figure Since the composition of the ashes from the H/He burning is still unknown, we adopt a simple 
model in which we assume that prior to electron capture the gas consists of protons, electrons and a single 
species of nucleus of mass Anip and charge Ze. The mass fraction of protons X is given by pp = Xp, so that 
the density of nuclei is pN = i^~X)p and the local number ratio of protons to nuclei is Up/ tin — XA/{1—X). 
The mean molecular weight per electron for this mixture is 

^^^=X + ^(1-X), (8) 

p fie A 

where rig is the electron number density. We assume that only one nuclear species is present at a particular 
depth. Given the starting values Xr and Ai, the value of X at a particular depth fixes the mass of the 
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nucleus A due to baryon conservation, A/{1 ~ X) = const = Ai/{1 — Xr). Thus A, in some sense, represents 
the average weight of the nuclei present, and is a continuous variable. The models we discuss throughout this 
paper are for the "neutron-rich" case in which we neglect f3 decays and keep Z fixed and let A grow without 
bound. We have also constructed models in which the nuclei follow the valley of stability, i.e. at a particular 
depth we choose the most stable Z corresponding to the local value of A. Clearly, the real composition of the 
layer lies somewhere between these two cases. However, we find that (except for the density discontinuity 
mode, see §4.1) the adiabatic seismology depends very weakly on the exact nuclear physics in the layer, and 
so we adopt the "neutron rich" models for simplicity. We have yet to fully explore the rich nuclear physics 
of this transition layer. 



2.3. The Structure of the Hydrogen Electron Capture Boundary Layer 

At the accretion rates of interest for the bright X-ray sources, the downward diffusive drift speed of a 
nucleus in the hydrogen-rich region is much less than the downward flow speed from accretion, v = —vz = 
—mz/p. The nuclei thus do not have time to settle out (gravitationally separate) from the hydrogen before 
the electron captures occur. For example (using equation 4.3 from Bildsten, Salpeter & Wasserman 1993), 
for temperatures in excess of 5 x 10^ K, there is no relative diffusion of protons and heavy ions provided 
M > 10~^''AfQyr~"'^. Thus there are no relative separations in the continuity equations (Brown & Bildsten 
1998), and we construct models of the hydrogen electron capture boundary layer by integrating the continuity 
equation for protons, 

diij ^ 

^ + V ■ (upv) = -UpR,,. (9) 

We rewrite this equation in terms of the column density, y = J p{z)dz = p/g, and the hydrogen mass fraction 
X, giving 

dX dX 

-l^+^^ = -^Rec, 10 
ot ay 

where we have used the mass continuity equation. In a steady-state, the electron captures are balanced by 
the accretion flow, so that 

dX , , 

m— = -XR,,. 11 
dy 



Since the temperature gradient across the transition layer is small (see §5.1), we first presume the layer is 
isothermal in our calculations. In this case, equation (|l|) gives a complete description of the transition layer. 
We use these isothermal models in §^ to illustrate the spectrum of ocean g-modes. We discuss realistic non- 
isothermal models of the whole ocean in §^, in which the heat equation must also be integrated. A discussion 
of the equation of state and the microphysics of the ocean can be found in BC95. The temperature of the 
ocean is Tg = (T/IO^ K) « 5-7. 

We now make a simple model of the hydrogen electron capture layer at zero temperature. In this case 
no electron capture occurs until Ep > Q. We neglect the electron rest mass, assume the boundary layer is 
thin, write the pressure as p — gy = neEp/A^ and use the capture rate of equation (^). We then integrate 
equation ( pT| ) analytically, obtaining X[Ef) = X{Q)exp{—B), where 

R,r , 27r 1 In 2 ^ ^ ^ ^ 



^ rh ^ 9{hc)^ mg ( ~~ -2 ) ( ^ ^) ( ) 
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If we measure the characteristic "thickness", AEp = Ep — Q, of the electron capture layer by the depth at 
which the hydrogen is 90% depleted, we find 

/ • \ 1/4 / \^/^ / O \ -5/4 / ff \ 1/4 

^^^^^''"^^^fej ( 2xlOi4cms-0 (l^) [Tk-s) • ^''^ 

Integrating the equation of hydrostatic equilibrium, dz = —dy/p gives the corresponding physical thickness 
of the layer, 

^ 826cm / 9 \ ( Q \ ( Jl_\ 

E \faEddj V2 X 1014 cm s-2yl ^1293 keVy Vl065 sy ' ^ ' 

where 

1 Z„ „> , (-3/4)!^ , / I 



E -A<'-*' + ^*nSj.' '''' 

and the subscript i denotes initial values, is the initial hydrogen mass fraction, (—3/4)! = 3.63, and 
we have assumed that Z is constant through the layer. We include the scalings with ft and Q in these 
formulae as they are equally applicable to electron capture layers associated with nuclei at much larger 
depths. Numerical integrations of equation (^l|) with the full T = equation of state for the electrons and 
the exact expression for the electron capture rate (equation [^) find good agreement with both the scalings 
and prefactors given in equations (|3|) and (<; 5% for m <^ lOrriEdd after which the term oc {Ep — QY 
neglected in equation (^ becomes important). 

As shown in Figure H, at high enough temperatures significant electron capture occurs pre-threshold. 
For a particular m, there is a critical temperature Tc above which this occurs, roughly given by ksTc AEp 
from equation (^3|). We accurately estimate Tc by integrating the continuity equation ( |ll|) with the pre- 
threshold expression for the electron capture rate (equation (H). Again neglecting the electron mass, we find 
that when 

^'^^^"^■n^y Ux 10i4cms-2j (,1293 kevj ' ^^^^ 

more than 50% of the protons have captured electrons pre-threshold. The prcfactor in this expression is 
taken from numerical integration of equation ( [III ) for an isothermal layer with the full electron capture rate 
(eq. 1^) and equation of state. Figure ^ shows the critical lines in the m-T plane for which 50% and 90% 
of the protons capture pre-threshold. There is excellent agreement with the analytic estimate. 



3. Adiabatic Oscillations in the Ocean 

Our present focus is on those g-modes that reside predominantly in the liquid ocean that lies beneath 
the H/He burning layer. The oscillations are adiabatic there since the local thermal time (hours to days) 
is much longer than the g-mode period. We first discuss the different sources of buoyancy in the ocean and 
then describe the adiabatic perturbation equations and our method for solving them. 



3.1. Sources of Buoyancy 



The g-mode frequencies are set by the local buoyancy force, which results from the density contrast 
between a displaced fluid element and the background fluid. There are two contributions to this density 
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contrast. The first is thermal buoyancy, which arises because the change in density of the fluid element as it 
is adiabatically displaced differs from the change in the background density due to the non-isentropic density 
gradient. The second contribution is from composition gradients. For example, in the hydrogen electron 
capture layer the timescale for electron capture is much greater than the mode period, so that the perturbed 
fluid element's composition is fixed during the pulsation. Because of the large composition gradient in the 
layer, the perturbed element is heavier than the surroundings and is forced back. 

The local buoyancy is measured by the Brunt- Vaisala frequency N, given by 

A fdhjp 1 d\np\ /d\np^ 1 \ ^^^^ 

\ dz Ti dz J ' \ dz ri/iy 

where Fi = {dlnp/dln p)s is the adiabatic index and A is the convective discriminant. For a mixture of 
electrons and ions, in which the pressure is a function of density p, mean molecular weight per electron /ig 
and mean molecular weight per ion /i^, we can rewrite A as 

^^dlnp/J^ 1\ xrdlnT x^^dlnpe Xt^.dlnpi 

dz \xp Ti/ Xp dz Xp dz Xp dz 

where xq = (91np/c)lnQ) with the other independent variables held constant. Using the thermodynamic 
identities Fi = (Fg - 1)xt + Xp, and V ad = {dhiT / d\np), = F3 - 1/Fi (Cox & Giuh 1968), we obtain 



N^h _ XT 
9 Xp 



Va, 



d 



d\nT 
dlnp 



X^e f d\np,e \ _X4M f d\np,i \ 
Xp \ dlnp J ^ Xp \ dlnp ) ^ 



(19) 



where the subscript * refers to the stellar model and h = p/pg is the local pressure scale height. The first 
term in this equation is the thermal buoyancy and the remaining terms are contributions from composition 
gradients. Figure ^ shows the different contributions to in isothermal models of the ocean. In the electron 
capture layer, the buoyancy due to the /ie gradient dominates the other contributions, elsewhere the thermal 
buoyancy is most important. The buoyancy due to the pi gradient is of order the thermal buoyancy in the 
layer; despite Xpi being small (the pressure is insensitive to /Zj) the change in pi is typically large across 
the layer, Api/pi w 6. We use isothermal models here for illustrative purposes. In reality, the ocean has a 
slight temperature gradient (see j |5.1| ), which changes the thermal buoyancy. However, the electron capture 
transition layer is typically thin enough so that it can be described by a single temperature. 

We now make an analytic estimate of the thermal buoyancy in the isothermal case under some simplifying 
assumptions. First, we use the relation Vad = xtp/ pTcyTi to write the thermal buoyancy as 

Nl^^Af;, (20) 

where cy is the heat capacity at constant volume. There are two contributions to xt- At low temper- 
atures, the electron contribution is small, so that xt is dominated by the contribution from the ions, 
XT pksT / piTUpP as T — > (see Hansen & Kawaler 1994 for a useful discussion). The heat capacity 
is 

Pimp \2 pe hp J 

where the first (second) term is the contribution from the ions (electrons). For this estimate, we ignore the 
electron contribution and treat the ions as an ideal gas, neglecting Coulomb corrections to the equation of 
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state. Our detailed numerical calculations include all of these effects. Taking Fi « Xp ~ 4/3 for relativistic 
electrons, the final result is (compare BC95 equation 3.9) 



Nik',. ^ (22) 

This estimate is shown as the dot-dashed line in Figure ^. It performs well at low temperatures, but under- 
estimates the thermal buoyancy at higher temperatures, where the contribution from the electron entropy 
becomes non-negligible. The non-ideal gas corrections for the ions (for which we use the parametrization of 
Farouki & Hamaguchi 1993) have little effect on N. 

At zero temperature, the thermal buoyancy vanishes and all buoyancy comes from the /ig gradient. 
The pressure is entirely from the degenerate electrons, giving Fi = Xp = — X/^e ^^'^ X/^i = 0, so that 
N'^ = —g din fie/dz. This form of buoyancy has been studied before in the core of a neutron star where 
the density dependent /3-equilibrium means that fXe changes with position. Reisenegger & Goldreich (1992) 
found a new class of core g-modes in the <^ 100 Hz range from this effect. These modes do not penetrate far 
into the crust and so most likely have little impact in the neutron star ocean which we are studying. 



3.2. Solving the Adiabatic Mode Equations 

The equations describing linear adiabatic perturbations of the thin (h <ti R) ocean are 

d^z Sp 
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(23) 
(24) 



dz p 

d 5p _5pl l \ , iz {(^'^ N'^^ 

dz p p h 

where is the radial displacement of the perturbed fluid element and dp is the Eulerian pressure pertur- 
bation. We follow BC95 and use the notation S to describe an Eulerian perturbation and A to describe a 
Lagrangian perturbation, A = 5 -I- ^ • V. We adopt the Cowling approximation, in which we neglect per- 
turbations of the local gravity. We use Newtonian gravity and take g = GM/R^. These equations are in 
the plane-parallel limit since the depth of the ocean is much less than the stellar radius. The horizontal 
wavenumber is then k = y/UJ + l)/ R for a slowly rotating star. Bildsten et al. (1996) showed that rapid 
rotation simply modifies this angular eigenvalue, allowing our calculations to be easily extended to the case 
of rapid rotation (/<, ^ mode frequency). 

Solving the adiabatic mode equations requires h, Fi and N throughout the ocean. BC95 showed that 
the shear modulus of the crust effectively excludes g-modes because they have large transverse shearing 
motions. We thus demand that the radial displacement vanishes, = 0, at the ocean floor where the ions 
crystallize. This is defined by F = F„i « 173 (Farouki & Hamaguchi 1993), where 



measures the importance of Coulomb effects. Here a is the average ion spacing, Aira^ni/^ = 1, and we 
presume for simplicity that only a single ion species is present. Assuming relativistic electrons and using the 
T = Fermi energy (very good approximations at the crystallization depth), we find that the crust starts 
at a column density 

... 1.8. X (f )' L^jy (§y"" (^t (26) 
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The ocean can crystallize before hydrogen electron captures occur only if Tg < 1 (Z/SO)^/"^, an unlikely 
circumstance at the accretion rates we are considering. 

We adopt a "shooting" method for solving the eigenvalue problem. For a particular trial frequency, 
we integrate from the surface of the ocean to the top of the crust. The boundary condition at the ocean 
floor, ^2 = 0, is only satisfied when the trial frequency is an eigenvalue of the equations. We start our 
integrations at an arbitrary place (column depth yt = 10^ g cm~^) beneath the H/He burning layer. We 
apply the boundary condition there that the Lagrangian pressure perturbation vanishes [Ap — 0). Although 
not strictly correct, this boundary condition is reasonable for the low frequency g-modes we discuss, and we 
find that the mode frequencies are insensitive to changes in the upper botindary condition or starting point 
for the intcgrations.n 



4. The Spectrum of Ocean g-Modes 

The distinct sources of buoyancy in the region beneath the H/He burning yields an extremely rich 
spectrum of g-modes. The abrupt rise in density from the hydrogen electron captures supports a single density 
discontinuity mode, and the internal buoyancy due to the composition gradient creates a new spectrum of 
modes which are "trapped" in the transition layer. There is also a set of thermal modes, as originally 
described by BC95, which we now find are separated by the density discontinuity and confined to either the 
upper or lower parts of the ocean. 

We use isothermal models with Ai ~ 60 in this section in order to illustrate and explain the physics of the 
different types of ocean g-modes. The large amount of overlapping frequencies for these modes means that 
we must first understand them separately before jumping in to the realistic models that include temperature 
gradients. Our description of the mode spectrum of non-isothermal models in ^ is then much simpler. 
These stars also have a surface wave or f-mode and acoustic waves or p-modes. These modes, which have 
frequencies ^ 10 kHz, penetrate into the crust because they have little transverse shear and arc thus not 
well represented by our thin ocean models and so we do not discuss them further. 



4.1. Discontinuity Modes 

As Finn (1987) noted in his original study of discontinuity modes in isolated neutron stars, we gain 
some intuition by first considering waves at the interface of two constant density, incompressible fluids. This 
case can be solved analytically. Since the discontinuity lies closer to the surface of the ocean than the floor, 
we take the upper fluid (density p+) to have a finite extent H and the lower fluid (density to extend to 
infinite depth. In the shallow water wave limit {kH <C 1) the discontinuity mode frequency is then 

cj^gen^^-^. (27) 

p- 



^ The g-mode frequencies are insensitive to details of the upper boundary, despite there being some mode amplitude there, 
because they are mainly determined by buoyancy deeper in the ocean. The exceptions are the upper thermal g-modes, and at 
high temperatures the trapped modes, whose frequencies depend on the thermal buoyancy above the electron capture layer. 
For example, changing yt from 10® g cm"'^ to 5 X 10* g cm~^ changes the discontinuity and lower thermal mode frequencies 
by 1%, the trapped mode frequencies by up to 10% depending on the temperature, and the upper thermal mode frequencies 
by 20%. We di scus s the effe ct on the upper thermal modes and trapped modes of extending our calculations into the upper 
atmosphere in §4.2 and § [l.3| . 
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We apply this by ignoring the finite thickness of the hydrogen electron capture layer for now (sec §4.1.1 for 
finite thickness models) , and considering a simple model, a T = ocean with a discontinuity in p and /ie at 
the place where Ep = Q. The buoyancy then vanishes everywhere, except at the discontinuity. To solve this 
problem, we modified our shooting scheme and applied the matching conditions = S^^- and Ap_|_ = Ap_ 
across the discontinuity (Finn 1987; McDermott 1990). In these equations, a + (— ) refers to a quantity in 
the upper (lower) part of the ocean. These conditions are physically motivated; must be continuous at the 
discontinuity because of the continuity equation, while Ap must be continuous to avoid infinite acceleration 
of a fluid element. There is no restriction on the transverse displacement and indeed we find a large 
jump in S^±_ across the discontinuity. 

Equation ( p7| ) motivates a frequency estimate of the density discontinuity mode associated with the 
electron capture layer. The length scale governing the behavior of the eigenfunction in the compressible case 
is the pressure scale height just above the discontinuity 

-1 
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(28) 



where we neglected the electron mass and assumed relativistic degeneracy. A natural guess for the disconti- 
nuity mode frequency is then 

(29) 
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and at fixed pressure {p/ ^,e)+ = [p/ p,e)-, so that (at T=0) 
/d = lllHz' ' ' ' 



1 

Me/ + 



10 km^| + 1) 
R 



1/2 



or for our case of a single ion species. 
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(31) 



where AZ and A^ are the changes in the charge and mass of the nuclei from one side of the discontinuity 
to the other. The prefactors in equations (^) and (|3^) are from our numerical calculations. The scalings 
found are as predicted by equation ( p9[ ) . The mode frequency is insensitive to the position of the surface or 
floor of the ocean as long as they lie more than a few scale heights away from the discontinuity. This is as we 
expect in the spirit of our guess above; it is the shortest length scale near the discontinuity which dominates 
the mode eigenfunctions and sets the mode frequency in the limit of kH <C 1. 

The left panel in Figure ^ displays the eigenfunctions for a T = density discontinuity mode with 
Xr = 0.1. The density profile used for this mode is the dotted line in Figure |. There is a large contrast (and 
actually a sign change) in the amplitude of the transverse displacement above and below the discontinuity. 
This means that almost all of the mode energy lies in the upper part of the ocean, above the boundary layer. 
The radial displacement has maximum amplitude at the discontinuity, and decays towards the ocean surface 
and floor. 



4- 1.1. Discontinuity Modes with a Finite- Thickness Transition Layer 



We now consider models with a finite thickness hydrogen electron capture layer that are appropriate for 
the accretion rates we are considering. In these models, the density continuously changes across the layer. 
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The solid lines in Figure |^ show the density as a function of position for three different accretion rates when 
Tg = 5 and Xr = 0.1. The dashed lines show the changing mass fraction of hydrogen. The dotted line 
shows the case where all the captures occur at Ep = Q and T — 0. Clearly, such a sharp discontinuity is no 
longer present at these high accretion rates. However, we find that the density discontinuity mode has the 
same frequency as the discontinuous model as long as the electron captures occur over less than a pressure 
scale height. The eigenfunction is also qualitatively the same, except that the boundary layer is resolved, 
for example there is a large transverse shear as rapidly decreases^, and the radial node which at T = 
lies immediately below the discontinuity is moved deeper in the ocean. The right panel of Figure ^ shows a 
typical discontinuity mode in an isothermal ocean. 

It might come as a surprise that the same discontinuity mode appears in both the simple model and 
the models with a finite thickness electron capture layer. Indeed, one case involves matching conditions 
across the discontinuity while the other involves smoothly integrating the adiabatic mode equations through 
the layer. We now show that, in the limit of a thin layer, the matching conditions can be recovered by 
integrating the mode equations through the electron capture layer. We use the transverse momentum 
equation, ikSp — puj'^S,± to transform the dependent variables in the adiabatic mode equations (|2^) and ( p^ ) 
from {Sp/p,^z) to (C_L,^z), finding 
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where Cs is the adiabatic sound speed, = Tigh. The discontinuity mode has lo^ <C iV^, 
^2 ^ fc^^-Li allowing equation d32) to be approximated as 



(32) 
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Integrating this through the transition layer gives 



l + O 



(35) 



where Az <^ /i is the thickness of the layer. The pressure is nearly constant within the electron capture layer 
as p and /ig trade against each other. In other words the density scale height is much less than the pressure 
scale height, h{dhip/dz) 3> 1, so we simplify 
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(36) 



and equation (B3f) becomes 



d£,i_ _ d\np / ^ ^ .kg ^ 
dz dz \ up- " 



(37) 



^ One might worry that the large transverse shear in the transition layer would be unstable to a Kelvin-Helmholtz type 
instability. At the accretion rates and Xr's we consider, however, we find that the shear in the transition layer is stabilized by 
the buoyancy there. The Richardson number in the layer is always much greater than 1/4 as long as the transverse displacement 
at the top of the ocean is less than a stellar radius. 
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Using the fact that ^2 is nearly constant in the layer, we integrate though the layer using the transverse 
momentum equation to find 



h-Ap- 

But since the pressure is constant, so is ph, giving us 



P4 



/ Az 
l + O — 



(39) 



For a thin boundary layer (Az <^ h), equations ( Pq ) and ( p9| ) are just the matching conditions from our 
simple model. 



4-. 1.2. Dependence of fd on Temperature and Accretion Rate 

The discontinuity mode frequency differs slightly from the simple estimate above (equation pl| ) due to 
three different effects: (1) pre-threshold electron captures change the position of the layer in the ocean, (2) 
thermal buoyancy provides some extra restoring force above and below the transition layer, and (3) at high 
TO and T the thickness of the layer can exceed the local scale height. However, we found that the dependence 
on these factors is rather weak and typically is not larger than 20-30% for an order of magnitude change in 
TO or T. 

Figure ^ shows the effect of to and T on fd for isothermal oceans with Xr = 0.1. At low accretion rates 
(to <^ 0.05TOEdd)j increasing Tg decreases fd because pre-threshold captures move the electron capture layer 
upwards in the ocean, decreasing the scale height at the transition layer. At higher m's, fd increases with Tg 
because of increasing thermal buoyancy and thickness of the transition layer. The effect of thermal buoyancy 
is shown by the dotted lines, which show fd for models where the thermal buoyancy is omitted. The effect of 
increasing layer thickness is shown by the T — line, which can be fit by a simple m-^^^ scaling, as suggested 
by equation (^^, for the physical width of the electron capture layer. In the relevant temperature range of 
Tg « 5-7, we find that the frequency scales roughly as fd oc m^^^^. 

The models shown in Figure H are equilibrium models of the ocean. In reality, the electron capture layer 
is unlikely to be in complete thermal balance with the current value of the hydrogen mass fraction, as the 
thermal time there (minutes to hours) is much less than the accretion time (hours to days). The response of 
the mode frequencies to changes in temperature on a short timescale will be determined by changes in the 
thermal buoyancy and not in the structure of the layer itself. 



4.2. Trapped Modes within the Transition Layer 

We have also found a new set of g-modes which reside within the electron capture layer. The restoring 
force for these modes comes mainly from the buoyancy within the layer due to the composition gradient, 
and so they exist even when T = 0. However, their existence depends on the finite thickness of the electron 
capture layer, which is a characteristic of the layers in accreting neutron stars. These modes were not found 
in the density discontinuity mode studies of isolated neutron stars (Finn 1987; McDermott 1990; Strohmayer 
1993), as in that environment the electron capture transitions were presumed (most likely safely) to be sharp. 
Though the modes propagate above and beyond the transition layer, we will refer to them as "trapped" 
modes, as we can understand their frequencies depending on how many nodes are present in the transition 
layer. 
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Figure ^ shows the ritr = 1 and ritr — 2 trapped modes in a T = ocean, where ntr is the number 
of radial nodes within the transition layer. Almost all of the mode energy lies in the transition layer in the 
unphysical T = limit. The transverse displacement is constant in the upper and lower oceans because N 
vanishes there and the mode is mostly transverse (see eq. |3^). When there are many radial nodes per scale 
height {kzh 3> 1) we can estimate the mode frequencies in the WKB limit. The coefficients in the adiabatic 
mode equations ( p3| ) and (|2^) are then constant and are integrated to obtain the dispersion relation for high 
order g-modes, — N'^k'^ /uP' . We require that there be an integer number of wavelengths, J kzdz = nir, 
giving the frequency of the n'th mode as 

fn~^- I Ndz^^- [ NhdXny. (40) 



27r2 n J 271^ n _ 

The trapped mode frequencies are set by the internal buoyancy due to the composition gradient in the 
transition layer. For nearly all cases we are considering, this dominates the thermal buoyancy within the 
layer (see Figure |5|) , so that N in that region is given by 

^2 g din lie g ^dlnX 
h d\ny h dlny ' 

where we have presumed (for simplicity) that the nuclear charge Z is constant so that fi^. is only a function 
of X. Adopting a relativistic equation of state, we use the solution X{Ef) = X{Q) exp(— i?), with B given 
by equation (|l^), to evaluate N . Integrating through the transition layer gives the convenient formula 

8.5Hz i^X^V^V rn g ft ^ 1293 keV 
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where we insert the prefactor found numerically. The analytic prefactor is close to that found numerically, 
and is a function of fundamental constants. We also show the scaling with ft and Q, since this formula is 
equally applicable to other electron capture boundary layers. 

We find excellent agreement between the predicted scalings and our numerical results. Figure ^ com- 
pares the analytic and numerical scalings with m and for a range of rit^ values. The left panel displays 
the first eight trapped mode frequencies as a function of X,. for a fixed accretion rate (rh = ra-E.dd) and T = 0. 
The filled triangles are our numerical results and the solid lines are from equation (^2|). The right panel of 
Figure |lO| compares the m scaling for Xr = 0.1 and T = 0. It is clear that the WKB formulation gives an 
excellent answer at high ntr and is reasonable at low ntr as well. We close by noting that the ratio of the 
discontinuity mode frequency to a trapped mode frequencies is nearly independent of everything 
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(43) 
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and thus might be an excellent observational discriminant. 

For T > 0, the effect of thermal buoyancy is to allow these modes to extend into the upper ocean and 
atmosphere. Figure ^ shows the ntr = 1 and ntr = 2 trapped modes in an isothermal ocean with Tg — 5. 
The frequencies are not much different from their T = values. The transverse displacement in the electron 
capture layer is less than at T = because the thermal buoyancy above the layer makes decay with depth. 
The mode extension into the upper ocean is critical for excitation of the trapped modes, particularly as the 
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mode then extends into the H/He burning layer. We are currently investigating this possibility, as well as 
excitation from the electron captures themselves. The thermal buoyancy is also larger at lower pressures 
and this will increase the mode frequency when the mode amplitude is significant in the upper atmosphere. 
We now show that at finite temperatures there is the additional complication due to the coincidence that 
the thermal modes are "mixed" with these trapped modes, at least in frequency space. 



4.3. Thermal Modes 

The thermal buoyancy throughout the deep ocean of degenerate electrons produces a set of thermal 
g-modes (BC95). These modes were extensively discussed by BC95, who pointed out the weak frequency 
dependence on the depth of the ocean, a convenient situation since the depth depends so strongly on uncertain 
quantities (see eq. [^). Just for completeness, we estimate the thermal g-mode frequencies in the WKB 
limit by integrating equation (po|). In a region where /i^ is constant and the electrons are very degenerate, 
equation ( p2|) gives N'^h'^ « diksT /i^iirip = constant and the integration simply gives 
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where yt {yb) is the column depth at the top (bottom) of the ocean. The frequencies depends most strongly 
on the temperature and ion mean molecular weight and only logarithmically on the physical thickness of the 



We find that the density discontinuity due to the hydrogen electron captures divides the thermal g-modes 
into two types, which we call the upper and lower thermal modes. Examples are shown in Figure |l^ for an 
isothermal ocean with Tg = 5, m = mEdd and — 0.1. The upper thermal modes (left panel in Figure [l^ ) 
have most of their energy above the electron capture layer. The layer acts as a floor for these modes and 
there is a node in the radial displacement there, just as there would be if we demanded our normal bottom 
boundary condition. This is not by design, but is rather how the modes behave. This motivates us to use 
our WKB estimate (eq. Q]) and insert typical values for jii and column depths in order to estimate these 
upper thermal mode frequencies 



_ 9.0Hz /Tg QV'\ 

/th.uppcr - It-—' 



Vb 



2 X 1010 g cm-2 J VlO^ g cm-2 



yt 



(45) 



10 km\ fl{l + 1) 



1/2 



R J \ 2 

The lower thermal modes (example in the right panel of Figure |l^) have most of their energy beneath the 
electron capture layer. We find that the Lagrangian pressure perturbation has a node at the electron capture 
layer, just as it would if we put our traditional top boundary condition at that location. In that sense, the 
layer acts as a surface for these modes. Equation (H) then gives the lower thermal mode frequencies as 
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where we have put in the top column as the electron capture layer and the bottom at the crystallization 
depth. We also took /i^ = 60, as no hydrogen is present at these depths. The upper and lower thermal mode 
frequencies are similar because the difference in /i^ above and below the electron capture layer is compensated 
by the greater depth of the ocean beneath the electron layer than above. Here niower (^upper) is the number 
of radial nodes below (above) the electron capture transition layer. 

The thermal mode spectrum is split because the jump in density and A'^^ across the transition layer 
creates an impedance mismatch, so that an impinging wave is partially reflected by the transition layer. This 
phenomenon is seen in so-called "mode trapping" in white dwarfs, in which thermal modes are "trapped" in 
the thin hydrogen or helium layers at the surface of the star (Winget, Van Horn & Hansen 1981; Brassard 
et al. 1992). These modes (analogous to our upper thermal modes) are believed to be preferentially excited. 
It may well be that something similar happens in our case. Just like the trapped modes (§4.2), the upper 
thermal modes have significant mode energy near the upper boundary, and so it is likely that they will 
extend upwards into the upper atmosphere. We are currently investigating this as we consider excitation 
mechanisms. The thermal buoyancy in the upper atmosphere will most likely increase the frequency of the 
upper thermal modes. For example, Strohmayer & Lee (1996) investigated the seismology of an atmosphere 
accreting and burning in steady state and found thermal g-modes of frequency approximately 50Hz for 
comparable accretion rates. The lower thermal modes are more or less confined to the deeper regions of the 
star. They do propagate upwards, but not enough to have their frequencies modified by any atmospheric 
physics. 



4.4. Avoided Crossings and Mode Identification 

To conclude our discussion of isothermal models, we show in Figure ^ the spectrum oi I = 1 modes in 
an isothermal ocean with Tg = 5, m = rriEdd and Xr = 0.1, our fiducial isothermal model. Each horizontal 
line shows a mode frequency. The highest frequency mode is the discontinuity mode, lower frequency modes 
are trapped and thermal modes. The positions of the radial nodes in each mode are shown as circles and 
triangles. The dashed vertical lines approximately bound the electron capture layer, showing where the 
thermal and composition gradient contributions to N'^ are equal. 

Diagrams like the one in Figure |l^ help greatly in mode classification. The type of a given mode can be 
inferred more or less by where the extra node is placed as one jumps from mode to mode down the diagram. 
Trapped modes have a new node added within the transition layer and this node remains there as one moves 
down the diagram. These nodes are marked with circles. A lower (upper) thermal mode has a new node 
placed at the bottom (top) of the transition layer and this node then moves downwards (upwards) in the 
ocean as one moves down the diagram. These nodes are marked with downward-pointing (upward-pointing) 
triangles. Thus in Figure the types of modes from top to bottom are d,u,t,l,u,l,t,l,u,t,l where d is for 
discontinuity, t is for trapped and u(l) is for an upper (lower) thermal mode. This method of classifying 
nodes greatly eases the mode identification problem. In particular, it helps to distinguish between trapped 
and upper thermal modes, whose eigenfunctions may be very similar^ (for example, compare Figures |ll|and 
p^ ). It is not completely reliable however, as it is complicated by the mixing of mode eigenfunctions near 
avoided crossings, as we now discuss. 



*This is because the thermal buoyancy above the transition layer begins to play a role in setting the trapped mode frequencies. 
The trapped and upper thermal modes remain two distinct sets of modes, however. For example, their frequencies depend 
differently on the position of the upper boundary. 
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It is by complete coincidence that the upper and lower thermal modes are at a similar frequency to the 
trapped modes. The frequencies of these modes depend differently on m, Tg: and X,., leading to "avoided 
crossings" whenever two mode frequencies try to cross and an accidental degeneracy arises. Avoided crossings 
were found for non-radial oscillations in massive main sequence stars by Osaki (1975), and in the neutron 
star context by Carroll et al. (1986), who found avoided crossings between g-modes and p-modes in their 
investigation of non-radial oscillations in the presence of strong magnetic fields {B ^ 10^^ G). 

Avoided crossings can be seen in Figure |lj which shows the frequencies of the discontinuity mode (upper 
solid line) and several thermal and trapped modes as a function of Xr for an isothermal ocean with Tg — 5 
and m = mEdd- Each solid line is the frequency of a mode with a fixed number of radial nodes n, and, as 
expected, the frequency always decreases as n increases. The different dependencies of the mode frequencies 
on Xr result in a series of avoided crossings. Along a solid line, the number of radial nodes is conserved, 
but the character of the mode (upper thermal, lower thermal or trapped, denoted by an upward-pointing 

1/2 

triangle, downward-pointing triangle or circle) may change. The dot-dashed line shows the Xr scaling 

— 1/2 

expected for the trapped and discontinuity mode frequencies. The thermal mode frequencies scale as /Xj . 
For the lower thermal modes (dashed line) this gives a decreasing frequency with Xr because fii in the lower 

1 /2 

ocean increases with Xr- For the upper thermal modes (dotted line) the scaling is oc Xr at high Xr but 
levels out at low Xr- On average, the different types of modes follow the predicted scalings by changing n 
in avoided crossings. Locally, however, the behaviour with Xr may be quite different. 

Figure ^ is a similar diagram to Figure |l^ in which we show the dependence of the mode frequencies 
on temperature in an isothermal ocean with m = mEdd and Xr = 0.1. Again, each solid line is for a mode 
with a fixed number of radial nodes n, and, as before, the different scalings of the mode frequencies with 
temperature yield avoided crossings. The dot-dashed line shows the T^/^ frequency scaling of the thermal 
modes. The dashed line shows the frequency of the n = 1 trapped mode when thermal buoyancy is not 
included. It shows how the frequency of the trapped mode changes due to changes in the structure of the 
boundary layer. Again, on average the modes follow the expected scaling with T, but not necessarily locally. 

Avoided crossings have been discussed by Aizenman et al. (1977), Gabriel (1980) and Christensen- 
Dalsgaard (1981) using an approach analogous to degenerate perturbation theory in quantum mechanics. 
In this picture the frequency splitting at the point of closest approach during a crossing is proportional to 
the "overlap" between the eigenfunctions of the two modes (the off-diagonal matrix element in quantum 
mechanics). This may explain why an avoided crossing between an upper thermal mode and a trapped mode 
has a larger frequency splitting than an avoided crossing involving a lower thermal mode in Figures ^ and 
|l5|. The "overlap" between upper thermal mode and trapped mode eigenfunctions is greater than between 
an upper thermal or trapped mode and a lower thermal mode. 

During an avoided crossing, the mode eigenfunctions are mixed. We show this in Figure ^ in which 
we take the avoided crossing which is circled in Figure |lj, and show the energy density of the two modes 
before, during and after the avoided crossing. This crossing is between a lower thermal mode and a trapped 
mode. The n = 4 and n — 5 mode energy densities are shown as a function of column depth for a range of 
temperatures from top to bottom. The circles show the position of the radial nodes. At Tg = 1.3, the lower 
frequency n = 5 mode is a lower thermal mode, while the higher frequency n = 4 mode is a trapped mode. 
At the avoided crossing, Tg w 1.7, the mode eigenfunctions are qualitatively similar, differing only by a node. 
At higher temperature, Tg = 2.3, the two modes have exchanged characters. The lower frequency n = 5 
mode is now a trapped mode while the higher frequency n = 4 mode is a lower thermal mode. The mixing 
of mode eigenfunctions may have implications for mode excitation, allowing exchange of energy between the 
two modes (for example, see Christensen-Dalsgaard 1981). 
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5. Non-Isothermal Ocean Models 

Up to this point we have presumed that the ocean is isothermal and considered the temperature to be a 
free parameter. We now construct non-isothermal models of the ocean by integrating the heat equation (see 
Brown & Bildsten 1998) through the ocean. There are several sources of energy that contribute to the heat 
flux in the ocean. The first is the energy released when the neutrons from the hydrogen electron captures 
combine with the heavy nuclei, roughly Eh w 7 x 10^^ ergs g^-'^ « 7 MeV per accreted nucleon. The second 
is the energy released as matter falls downwards (or what is sometimes called "compressional heating"), 
-Ecomp ^ ksT / ^iVTLp. Finally, there is energy released in the deep crust by nuclear electron captures and 
pycnonuclear reactions. Most of this energy (about an MeV per accreted nucleon, Haensel & Zdunik, 1990a) 
goes into the NS core and is lost via neutrinos, about 10% {Ecr ^ 0.1 MeV per accreted nucleon) flows 
upwards through the ocean (Brown & Bildsten 1998). The total heat flux above the electron captures (at 
say y k, l(f g cm~^) is 

Fa ccp ~ Th[Ef{J^j^ + -£^comp ^ -Fcr). (^7) 

The largest uncertainty in Jdccp is the amount of hydrogen left unburned and advected downward, Xr (Taam 
et al. 1996; Schatz et al. 1997b). As we show in this section (in agreement with Taam et al. 1996), for 
typical values of Xr ~ 0.1 and high accretion rates the region between the H/He burning layer and the 
hydrogen electron captures has a substantial temperature gradient. The other contributions to -Fdecp are 
typically not large enough to require a substantial temperature gradient (BC95) at sub-Eddington accretion 
rates, especially in the relativistic ocean underneath the hydrogen electron captures. 

Even when Jdccp is known and nuclear reactions beyond electron capture are ignored^, there are still 
additional uncertainties. One is the radiative opacity for this complicated mixture of elements. For pure 
iron at the depth where the hydrogen electron captures are occurring, the predominant heat transport 
mechanism would be electron conduction (Gudmundsson, Pethick, & Epstein 1983). Thus for simplicity, 
we neglect radiative heat transport and assume that all the heat is carried by electron conduction. For the 
conductive opacities, we use the results of Yakovlev & Urpin (1980) for electron-ion collisions, and the fit of 
Potekhin et al. (1997) for electron-electron collisions. An additional uncertainty is the choice of temperature 
at the top of the ocean when the H/He burning is time dependent. This is because the separation between 
the electron capture region and the burning location is not far enough that a simple radiative-zero like 
solution applies, at least when T > 5 x lO^K at the electron capture depth. 



5.1. A Few Illustrative Models 

In order to illustrate the effects of -Fkccp, we show in Figure ^ the temperature, hydrogen mass fraction 
and density as a function of column depth for models with m = O.lmEdd and Xr = 0.2,0.1 and 0.05 (from 
top to bottom in the T and X panels). In the density panel, the curve with the highest density at the 
largest column depth has Xr = 0.2. These models have the outer boundary condition T = 2 x 10^ K at 
y = 10^ g cm~^ and, for simplicity, we set iScomp = and E^r = lO^^erg g^^ w 0.1 MeV per accreted 
nucleon. The initial nucleus has Ai = 60, Zi = 30 and is allowed to become arbitrarily neutron rich. The 



^For the purposes of illustration in this section, we persist with our simple model of electron capture followed by neutron 
capture onto a single type of nucleus. There is no reason to undertake a more complete treatment of the neutron captures//3- 
decays and the possibility of direct proton captures onto the nuclei when T > 10^ K until we have better knowledge of the 
nuclear mix from the H/He burning. 
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flux at depths far below the electron capture region is i^bottom ~ niEcr = 7.5 x 10^° erg s"'^ cm~^. It is clear 
that, for these values of Xr, the ocean is substantially hotter than if there were no hydrogen present. At this 
accretion rate, the ocean is nearly isothermal below the hydrogen electron captures. Most of the electron 
captures are post-threshold {Ep > Q) and occur 15-30 days after the matter has arrived on the star. 

Figure |8| shows models with a higher accretion rate, m — 0.5mEdd, and an outer boundary condition 
of T = 4 X 10^ K at y = 10^ g cm^^. The other variables are the same as in the previous example. The 
enhanced flux due to the higher accretion rate (see eq. |^) makes these models much hotter than the 
lower accretion rate ones. We may have slightly overestimated the temperature at these high m's, since for 
T > 10^ K, radiation may start to carry a significant fraction of the flux. Many of the electron captures 
for the Xr = 0.2 and Xr = 0.1 cases occur pre-threshold {Ep < Q). For the Xr = 0.2 case, one half of the 
hydrogen is depleted only 2.5 days after arriving on the NS. 

At these high temperatures (> 10® K) it is important to consider neutrino cooling. Using the formulae 
given by Schinder et al. (1987), we find that, even for the hottest model we consider here {m = 0.5, 
Xr — 0.2), neutrino cooling is unimportant at the depth of the hydrogen electron captures. It is certainly 
very important at greater depths however (Brown & Bildsten 1998). Once y ^ Ecrtn/ei,, where ei, is the 
neutrino cooling rate in erg g~^ s~^, the temperature gradient changes direction and heat starts to flow 
into the core, a possibility we neglect for now. For m = 0.5, this depth is y ~ lO^'^-lO^'* g cm~^, near the 
bottom of the ocean. At higher m's or X^'s neutrino cooling will become important at even lower column 
depths than this, although we stress again that at higher temperatures it becomes important to include the 
radiative opacity. 

How does the hydrogen electron capture transition layer compare to the isothermal cases we calculated 
previously? Figures |l^ and |l8| show that the transition layer itself is nearly isothermal. For example, even in 
the model with m — 0.5 and Xr = 0.2, the temperature changes by a factor of only 10% across the transition 
layer. We thus characterise the structure of the layer by the temperature at, say, the place where one half 
of the hydrogen is depleted. For example, the hydrogen mass fraction X in the non-isothermal model with 
TO = 0.5 and Xr = 0.1 agrees to within 10% at each depth with an isothermal model with T = 10® K and 
the same m and Xr- 

These examples are intended to illustrate the large impact of the hydrogen electron captures (Taam et 
al. 1996) on the deep structure of the neutron star. Much more work is needed on the products of H/He 
burning before any definitive results can be stated. However, we hope we have strengthened the case of 
Taam et al. (1996) that even small amounts of residual hydrogen at high m's can play an important role. 
We now describe the g-mode spectra of these non-isothermal models. 



5.2. g-Mode Spectra for Non-Isothermal Models 

The g-mode spectra of non-isothermal models of the ocean are easy to understand in the context 
of our study of isothermal ocean g-modes in The g-mode spectra for the to = 0.1, AT^ = 0.05 and 
TO = 0.1, ATr = 0.1 models discussed above are shown in Figure |l^. The discontinuity mode frequency is 
insensitive to details of the structure of the transition layer and depends mainly on Xr- It is given accurately 



by equation (31). The lower thermal mode frequencies are set by the thermal buoyancy in the ocean below 
the transition layer, which is almost isothermal. The frequencies are almost the same as the lower thermal 
mode frequencies in an isothermal ocean with a temperature equal to that just below the transition layer. 
The discontinuity and lower thermal mode frequencies are insensitive to details of the upper boundary and 
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thermal buoyancy above the transition layer. This is not true of the upper thermal and (to a lesser extent) 
the trapped modes. Hence we expect that their frequencies will change as we extend our models into the 
upper atmosphere. Also, we have not included radiative opacity which becomes important near the upper 
boundary and will reduce the temperature gradient there. In Figure |l^ we show the predicted frequency 
increases due to temperature for the thermal modes (T^^^ scaling) and due to X,. for the trapped modes 
(X^/^) by vertical bars. These simple scalings predict the difference in frequencies between the left and right 
panels of Figure |l^ fairly well. 

As a further example, we have reconstructed one of the models of Taam et al. (1996). They followed the 
time dependent evolution of models with accretion rates ranging from 0.1-1 times the Eddington rate. All 
of these models had substantial residual hydrogen, ranging at the highest m's from Xr = 0.08-0.21. Their 
neutron star had M = 1.4Mq and R = 9.1 km and the initial heavy nucleus was iron {Ai = 56, Zi = 26). We 
chose model 8 from their tables 1 & 2 which has Xr = 0.137 and m = O.lmEddj and varied our outer boundary 
condition until we matched their stated deep temperature. We then calculated the g-mode spectrum of this 
model. This is shown in Figure ^ 



6. Conclusions 

We have exhaustively studied how the hydrogen electron capture layer expected in the rapidly accreting 
{M > 10~^'^Mq yr~^) neutron stars affects the internal g-mode spectrum of these objects. The abrupt rise 



in density across the layer supports a density discontinuity mode (see §4.1) of frequency 



where AZ and are the average change in the nuclei's charge and mass through the layer due to neutron 
captures and subsequent /3-decays. The internal buoyancy due to the composition gradient within the 
electron capture boundary layer creates a new spectrum of "trapped" modes (see § |4.2| ) of frequency 

where ntr is the number of nodes in the boundary layer and we have omitted the weak dependence on the 
accretion rate. There is also a set of thermal g-modes in the same frequency range as the trapped modes, 
which are separated by the density discontinuity into two distinct sets (see §^). 

This work might eventually provide a natural explanation for sonic of the observed QPO's in accreting 
NS. Indeed, if the neutron stars in the "Z" sources were slowly rotating, the discontinuity (trapped) modes 
from the hydrogen electron capture layer would nicely overlap the freqencies of the QPOs seen in the 
horizontal (normal) branch. The stable 40 Hz QPO in the « 2 Hz X-ray pulsar GRO J1744-28 (Zhang 
et al. 1996) is an intriguing one to consider, as the discontinuity mode is in the right frequency range and 
might have a small enough shear so that the < 10^^ G field of this pulsar (Bildsten & Brown 1997) will not 
affect the mode frequency. However, there is still much theoretical work to be done, from understanding 
how the modes are excited to how they modulate the X-ray flux. 
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6.1. Observational Tests 

Observational progress will come by identifying QPOs with modes of different Z's, where the ratio of the 
mode frequencies are known. This is easiest when the neutron star is rotating slowly (spin frequency fs < 10 
Hz), as one then merely looks for the oc l{l + 1) scaling between QPO's that are (ideally) simultaneously 
present. The rapidly rotating case is more complicated as the angular eigenfunctions and dispersion relations 
are very different. As we discuss below, seismological progress can most likely be made only for those NS 
where the spin frequency has been measured during a Type I burst (see Bildsten 1998b for a summary of 
those objects). 

An important question is how fast can the mode frequencies change? The mode frequencies depend 
on internal conditions in the ocean, in particular on the temperature and X^- If a QPO was observed to 
change its frequency faster than the internal conditions in the ocean can change or if the dynamic range 
of the frequency was too great, this could rule out non-radial oscillations as a source of the periodicity. 
The discontinuity mode frequency and the trapped mode frequencies are sensitive to the residual hydrogen 
mass fraction, Xr- This changes from one X-ray burst to another, which will result in layers of different X^ 
being compressed towards the electron capture boundary layer. These layers are most likely Rayleigh- Taylor 
unstable and will mix, making the timescale for a large frequency change roughly the accretion time at the 
transition layer, days to weeks. The thermal modes and trapped modes are sensitive to the temperature in 
the ocean. The important timescale is the thermal time at the place where most of the mode energy resides, 
hours to days at the electron capture depth. Thus the upper thermal mode frequencies will change faster 
than the lower thermal mode frequencies. As we noted in since the thermal time is shorter than the 
accretion time, the change in the mode frequencies on a short timescale will be determined by the change in 
the thermal buoyancy and not the structure of the transition layer. 



6.2. The Rapidly Rotating Case 

There is good reason to suspect that these neutron stars are rapidly rotating, as the prolonged accu- 
mulation of material at this rate will most likely spin up the star. The coherent periodicities during Type I 
X-ray bursts seem to indicate 250-500 Hz spin frequencies, which are much greater than the g-mode frequen- 
cies, but still small compared to the breakup frequency Qb ~ {GM / R^Y^^ (« 2 kHz for a 1AMq,R = 10 
km star). When f2 ^ i7f,, the unperturbed star is spherical and the centrifugal force can be neglected, in 
which case the primary difference in the momentum equations is the Coriolis force. As a result, the g-mode 
frequencies depart significantly from the cj^ oc /(/ + 1) scaling (Papaloizou & Pringle 1978). 

BUC96 made progress on this problem within what is called the "traditional approximation" , where 
the radial and transverse momentum equations separate and the resulting angular equation must be solved 
to find the angular eigenfunctions (no longer Yims) and transverse eigenvalues A = kR^ (no longer l{l + 1)). 
The radial equations are identical to the non-rotating case, so that if ujq is the eigenfrequency for the / = 1 
mode of a non-rotating star, then the oscillation frequency (in the rotating frame) at arbitrary spin is 
uj — ajo(A/2)^/^. These are then transferred into the observer's inertial frame via uji = uj ~ mfl, so that 
going from our non-rotating results to the rotating star is straightforward. However, nearly any frequency 
can be predicted in the absence of prior knowledge of the spin frequency, as the retrograde modes cover a 
large frequency range as the NS spin is varied. 

Robust predictions are, however, possible for those NSs where the spin frequency is known from a Type 
I burst. By far the best case to date is 4U 1728-34 (Strohmayer et al. 1996; Strohmayer, Zhang & Swank 
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1997) which has fg = 363 Hz. For a known spin frequency, BUC96 and Bildsten et. al. (1998) have 
shown that one can predict the observed frequencies as a function of the internal conditions (e.g. T or Xr). 
Bildsten et al. (1998) asked how the observed frequencies for different angular eigenmodes would change as 
the internal conditions change. They found that the prograde modes are all at relatively high frequencies, 
near a kHz. However, these mode frequencies did not exhibit the large dynamic range of the observed kHz 
QPOs. The splitting for some of the prograde modes was found to be nearly constant, but not equal to the 
spin frequency. Other prograde modes had varying differences and could potentially be applicable to those 
kHz QPO's that do not have constant frequency separation (Sco X-1, van der Klis et al. 1997, 4U 1608-52, 
Mendez et al. 1998). Of particular interest were the retrograde modes that appeared at / < 100 Hz. These 
showed a large dynamic range for only a small change in internal conditions. Bildsten et al. (1998) thus 
pointed to the possibility of explaining some of the QPOs seen in the < 100 Hz range from 4U 1728-34 and 
other Atoll sources. 

Bildsten et al. (1998) pointed out that confirming non-radial pulsations in a star with a previously 

measured spin frequency can come by identifying a measured low-frequency QPO with a particular value of 
the non-rotating frequency and then searching in the data for the higher frequency prograde modes. Such 
an exercise has yet to be carried out and could potentially bear fruit, as few Fourier bins would need to be 
searched. The dispersion relations are complicated enough that such exercises can only be done in concert 
with theory. At present, theoretical predictions can most likely be done at the 10% level. Further theoretical 
work to include general relativistic effects might be needed for the very rapidly rotating stars (/« > 500 Hz) 
before such an exercise could be carried out in detail, say at the 1% level. 

6.3. Future Work on the Capture Layer, Mode Excitation and Coupling to Accretion 

There is still much work needed on the physics of the hydrogen electron capture layer. For example, a 
more realistic calculation would follow several species of heavy nuclei and incorporate neutron capture cross 
sections and /3 decay half-lives for each. It is important to carry out some of these calculations, as there 
is a possible large-scale instability that might develop underneath the layer. After all the electron captures 
have occurred, the nuclei will, on average, be neutron rich. These nuclei advect downwards and eventually 
/3-decay (they typically are not Fermi blocked) on timescales much longer than days. This may cause a 
density inversion as the released electrons will decrease the density, possibly leading to a Rayleigh- Taylor 
instability underneath the hydrogen electron capture zone. 

Which modes does the star like to excite? We are presently working on the internal excitation and 
damping of the adiabatic modes we have presented here. Following the work of SL96, we are starting by 
considering the excitation of the upper thermal modes, trapped modes and discontinuity modes by the H/He 
burning layers. Our new understanding of the mode structure will modify the competition between excitation 
in the burning layers and damping elsewhere in the star that they found. We hope to show that a sub-set of 
the infinite spectrum of oscillations we have found here will be excited. Finding the steady-state amplitude 
of the oscillation is beyond the scope of these calculations since it would require going beyond linear order. 

How can a mode modulate the X-ray luminosity? There are a few ways one can imagine. One pointed 
out by Bildsten et al. (1998) for the rapidly rotating case is from direct interaction with the accretion disk, 
as the g-modes are confined to a narrow equatorial band that would lie in the same plane as the accretion 
disk. If the NS does not lie inside the last stable orbit, then the disk should run directly into the surface, 
allowing for potential periodic modulation of the accretion environment. An alternative is to invoke a weak 
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magnetic field. The mode frequencies are not modified when B < 10^ G (BC95), but the thermal modes 
have tremendous shear and will thus, at high altitudes, change their nature (Carroll et al. 1986). Whether 
the resulting outgoing waves can affect the accretion flow remains an open question. 

We thank Hendrik Schatz and Michael Wiescher for conversations regarding the nuclear physics, A. 
Muslimov for suggesting we check for shear instability, and Greg Ushomirsky and Ed Brown for many in- 
sightful conversations and comments on the manuscript. This work was supported by the NASA Astrophysics 
Theory Program through grant NAG 5-2819. L. B. was also supported by the Alfred P. Sloan Foundation. 
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Fig. 1. — A schematic of the atmosphere, ocean and crust of a neutron star accreting at a rapid rate. The 
accreted H/He burns at a column depth y « 10^ g cm~^ or a density p « 10^~^ g cm~^. The ashes from the 
burning form a degenerate liquid ocean of unburned hydrogen and heavy nuclei. Hydrogen electron captures 
occur once Ep = 1293 keV, or y w 10^*^ g crn^^, p « 10^ g cm^''. The thickness of the electron capture 
transition layer is roughly a scale height. The ocean crystallizes and forms the crust at y 10^^ g cm~^, 
p » 10^~^° g cm~^. The ocean is nearly isothermal beneath the electron capture transition layer with a 
temperature of 5-7 x 10* K. For purposes of comparison, the ocean thickness in relation to the radius is 
about the same as that on the Earth. 
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Fig. 2. — The comparison of the exact integral to the analytic fit for the hydrogen electron capture rates. 
The solid lines show the ratio of the exact integral over the electron energy (eq. |^]) to that determined by 
our simple analytic fit (eq. [^) for temperatures of Tg = 1,2,5,8 and 10. These curves are tabulated and 
later used to correct the simple formula. The vertical dashed line denotes where Ep — Q. 



-30- 




1000 1200 1400 1600 1800 

Electron Fermi Energy (keV) 



Fig. 3. — The proton lifetime to electron capture as a function of the electron Fermi energy. The solid lines 
show the lifetime of a proton (l/_Rec) as a function of the electron Fermi energy for a range of temperatures. 
Starting from the uppermost solid curve, we display the lifetime for Ts = 2, 3, 4. 5, 6. 7, 8 and 9. The vertical 
dashed line denotes where Ep = Q- The dotted curves are the rough measure of the time it takes matter to 
reach that depth for m = mEdd (lower dotted) and m = O.lrhEdd (upper dotted). 
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Fig. 4. — Pre-threshold vs. post-threshold electron capture on hydrogen. The m-T plane is divided into 
two regions: at high T or low m, most captures happen pre-threshold {Ep < Q); at low T or high to, most 
captures happen post-threshold {Ep > Q). The filled (open) triangles show points at which 50% (90%) of 
protons have captured pre-threshold. The solid hues show the T oc rv}/'^ scaling (equation 10). 
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Fig. 5. — The Brunt- Vaisala frequency for isothermal models accreting at to = TOEdd- The solid line shows 
N, and the different lines show individual contributions. The thermal buoyancy {dotted line) dominates 
above and below the electron capture boundary, while the buoyancy due to the /ie gradient {short dashed 
line) dominates in the boundary layer. The Hi buoyancy {long dashed line) is about the same as the thermal 
buoyancy in the layer. The dot-dashed line shows the analytic estimate of the thermal buoyancy (equation 
(E3)), which performs best at low temperatures. 
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Fig. 6. — Eigenfunctions of the I = 1 density discontinuity modes. The left panel shows the density 
discontinuity mode for a T = ocean with Xr = 0.1. The right panel is for an isothermal ocean with 

rh = rrLEdd, Ts = 5 and Xr = 0.1. The middle (bottom) panels show the transverse (radial) displacement as 
a function of column depth. The top panels show the mode energy density (ii?32 /d logj^Q y {solid line) and 
Brunt Vaisala frequency AT on a linear scale {dashed line). The displacements shown are absolute values, so 
zero-crossings appear as cusps. We begin our integrations at a column depth y = 10^ g cm~^ and integrate to 
the floor of the ocean. In all our plots, we arbitrarily normalize the transverse displacement to = 10^ cm 
at the top of the ocean. The mode energy, Emodo is measured in units of E^2 = ^'mode/lO^^ erg. 
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Fig. 7. — Density as a function of depth in an isothermal ocean (Tg = 5) with an initial hydrogen mass 
fraction Xr — 0.1. The solid lines show (left to right) rn/ifiEdd — 0.01,0.1 and 1. The dashed lines (left to 
right) show the hydrogen mass fraction X as a function of depth on a linear scale, with the value at the left 
being Xr = 0.1. The dotted line is for a T = completely discontinuous model where all captures occur at 
Ef = Q- The density discontinuity mode for this density profile is shown in the left panel of Figure |[ 
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Fig. 8. — The effect of accretion rate and temperature on the I = 1 discontinuity mode frequency when 
Xr = 0.1. The soUd hnes show fd as a function of to for isothermal oceans with Tg = 0,1,5 and 10. 
The dotted Hnes show fd for models where the thermal buoyancy is omitted. At low to, fd decreases with 
increasing Tg. This is because pre-threshold captures move the electron capture layer upwards in the ocean, 
decreasing the scale height at the discontinuity. At high m, fa increases with Tg because of increasing 
thermal buoyancy and layer thickness. 
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Fig. 9. — The ritr = 1 and ritr = 2,1 = 1, trapped modes in a T = ocean with m = niEdd and Xr = 0.1. 
The restoring force for this mode is provided by the /ie gradient in the electron capture boundary layer. 
Almost all of the mode energy lies in the boundary layer. Labels and axes are the same as described in 
Figure ||. 
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Fig. 10. — The scalings of the I = 1 trapped modes with Xr and rh. The left panel shows the frequencies 
of the first eight trapped modes {ntr = 1 has the highest frequency) in an ocean with T = and m = niEdd 
{triangles) compared with the predicted scaling ftr = 8.5 Hz {Xr /O.lY^"^ /ntr {solid lines). The right panel 
compares the observed and predicted scaling with m, ftr = 8.5 Hz {rh/niEdd)^^^ /ntr, when Xr = 0.1 and 
T = 0. There is good agreement, even for small ritr where the WKB approximation doesn't strictly apply. 
Here ritr is the number of radial nodes within the transition layer. 
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Fig. 11. — The ntr — 1 and ntr ^ 2, I = I trapped modes in an ocean with Ts = 5, m = mEdd and Xr = 0.1. 
The effect of thermal buoyancy is to allow some of the mode energy to extend into the upper ocean. Labels 
and axes are the same as described in Figure |^. 




Fig. 12. — Thermal g-mode I — 1 eigenfunctions for an isothermal ocean with Ts = 5, m = mEdd and 
Xr — 0.1. Left (right) panel displays an upper (lower) thermal mode, so called as most of the mode energy 
resides above (below) the electron capture transition layer. Labels and axes are the same as described in 
Figure ||. 
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Fig. 13. — The spectrum of / = 1 modes in an isothermal ocean with Tg = 5, m = rfiEdd and Xr = 0.1. 
Each horizontal line shows a mode frequency. The highest frequency mode is the discontinuity mode, lower 
frequency modes are trapped and thermal modes. The dashed vertical lines show the places where the thermal 
and composition gradient contributions to are equal; these lines show the approximate boundaries of the 
electron capture layer. Moving down the diagram, each new radial node is marked according to the type 
of mode, circles for trapped modes and upward-pointing (downward-pointing) triangles for upper (lower) 
thermal modes. The trapped nodes remain within the boundary layer, whereas the nodes for the upper 
(lower) thermal modes move upwards (downwards) in the ocean as the frequency decreases. From the top 
down, the type of modes are d,u,t,I,u,l,t,I,u,t,I, where d is for discontinuity, t is for trapped and u(l) is for 
an upper (lower) thermal mode. See text for further discussion. 
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Fig. 14. — The I = 1 frequencies of the discontinuity mode and the first few thermal/trapped modes as 

a function of hydrogen mass fraction in an isothermal ocean with Tg = 5 and to = iriEdd- Each sohd 
hne is for a mode with a fixed number of radial nodes. The type of mode is shown as a circle (trapped 
modes) or triangle (thermal modes). The dot-dashed line is the x}^'^ scaling expected for the trapped and 

— 1/2 

discontinuity mode frequencies. The thermal mode frequencies scale as . For the lower thermal modes 

(dashed line) this gives a decreasing frequency with X^, while for the upper thermal modes (dotted line) the 

1 /2 

scaling is Xr at high X^ but levels out at low Xr- 
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Fig. 15. — The I — 1 frequencies of the first few thermal/trapped modes as a function of temperature in an 
isothermal ocean with to ~ rfiEdd and Xr = 0.1. As in Figure |lj, each sohd hne is for a mode with a fixed 
number of radial nodes, and we show the type of mode by a circle or triangle. The dot-dashed line shows 
the T^/^ frequency scaling of the thermal modes. The dashed line shows the frequency of the n — 1 trapped 
mode when thermal buoyancy is not included. The circled avoided crossing is displayed in detail in Figure 

m 
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Fig. 16. — The avoided crossing (for I = 1) between a trapped and lower thermal mode. The n — A and 
n — 5 mode energy densities are shown as a function of column depth for a range of temperatures from top 
to bottom. The circles show the position of the radial nodes. See text for discussion. This avoided crossing 
is circled in Figure Esl 
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Fig. 17. — The structure of an i? = 10 km, M = IAMq NS near the hydrogen electron capture region for 

an accretion rate of m = O.lmEdd and different values of Xr- The three hnes in the top (middle) panel show 
the temperature in units of 10* K (hydrogen mass fraction, X) as a function of the column depth y = p/g for 
(from top to bottom) Xr = 0.2, 0.1 and 0.05. The bottom panel displays the density (in units of 10*" g cm^'') 
for the three different values of Xr- The line with the highest density at the highest column is for = 0.2. 
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Fig. 18. — The NS structure near the hydrogen electron capture region for to = 0.5mEdd and different 
values of Xr- The labels and arrangement are the same as in Figure 
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Fig. 19. — The spectrum of ^ = 1 modes in non-isothermal oceans with to = O.lmEdd and Xr = 0.05 (left 
panel) and Xr — 0.1 (right panel). The temperature at the top {y = 10^ g cm^^) is Tg = 2 and rises to 
Tg « 4 {Xr = 0.05) and Tg w 5 {Xr = 0.1) below the electron capture layer (at y = 3 x 10^° g cm'^). 
We have included a deep flux of 0.1 MeV per accreted nucleon. The initial nucleus before electron capture 
has Ai — 60, Zi = 30. In each case, we show the discontinuity mode and the first ten trapped and thermal 



modes. The symbols in this Figure are the same as in Figure 13, For clarity, we have omitted the first 



radial node which lies near the upper boundary. The vertical bars show the difference in frequency expected 
between the left and right panels given the simple scalings and T^/^. 
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Fig. 20. — The spectrum of Z = 1 modes in a non- isothermal ocean with Xr = 0.137 and m = O.lmEdd- 
The temperature at the top (y — 10^ g cm"^) is Ts = 3.5 and rises to Ts = 5.6 below the electron capture 
layer (at y = 3 x 10^" g cm~^). We have included a deep flux of 0.1 MeV per accreted nucleon. The initial 
nucleus before electron capture is iron {Ai = 56, Zi = 26). This model is based on model 8 of Taam et al. 
(1996). We show the discontinuity mode and the first ten trapped and thermal modes. The symbols in this 
Figure are the same as in Figure |l^. For clarity, we have omitted the first radial node which lies near the 
upper boundary. From top to bottom, the modes are d,t,l,t,l,u,l,t,l,u,l, where d is for discontinuity, t is for 
trapped and u(l) is for an upper (lower) thermal mode. 



